{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这一节开始，介绍一些常用的概率分布（其实主要为了后续的LDA（隐狄利克雷分配）模型做准备），首先让我们从一枚硬币开始吧"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 一.伯努利分布\n",
    "假设我们手头有一枚硬币，用随机变量$x$描述扔硬币的结果，记$x=1$表示正面，$x=0$表示反面，另外$x=1$发生的概率记为$\\mu$，那么$x=0$的概率就为$1-\\mu$（$\\mu$未必为0.5），那么伯努利分布就是抛一次硬币的概率分布：    \n",
    "\n",
    "$$\n",
    "Bern(x\\mid \\mu)=\\mu^x(1-\\mu)^{1-x}\n",
    "$$   \n",
    "\n",
    "#### 均值、方差\n",
    "\n",
    "接着，我们可以很容易的求出该分布的均值、方差：    \n",
    "\n",
    "$$\n",
    "E[x]=1*u+0*(1-u)=u\\\\\n",
    "var[x]=E[(x-E[x])^2]=E[x^2]-(E[x])^2=1^2*u+0^2*(1-u)-u^2=u-u^2\n",
    "$$   \n",
    "\n",
    "#### 极大似然估计\n",
    "\n",
    "接下来，我们需要考虑考虑一个问题，假如给了我们一堆iid采样的样本$X=\\{x_1,x_2,...,x_N\\}$，如何去估计伯努利分布的参数$\\mu$，常用的一种方式就是极大似然估计，首先写出它的似然函数：    \n",
    "\n",
    "$$\n",
    "p(X\\mid\\mu)=\\prod_{i=1}^N[\\mu^{x_i}(1-\\mu)^{1-x_i}]\n",
    "$$   \n",
    "\n",
    "那么，其对数似然函数为：   \n",
    "\n",
    "$$\n",
    "ln[p(X\\mid\\mu)]=\\sum_{i=1}^N[x_iln\\mu+(1-x_i)ln(1-\\mu)]\n",
    "$$   \n",
    "\n",
    "令$ln[p(X\\mid\\mu)]$关于$\\mu$的导数为0，可求得：   \n",
    "\n",
    "$$\n",
    "\\mu_{ML}=\\frac{m}{N}\n",
    "$$   \n",
    "\n",
    "其中，$m=\\sum_{i=1}^Nx_i$。因此在极大似然估计的框架下，正面向上发生的概率是数据集里正面向上的观测所占的比例，但同时我们也可以发现一个问题，那就是如果观测量很少，极大似然估计的结果会出现**过拟合**的情况，比如我们连续抛一枚硬币3次，碰巧3次都是正面朝上，那么$N=m=3$，则$\\mu_{ML}=1$，这显然不太合理，后面我们会引入关于$\\mu$的先验概率来得到一个更加合理的结果"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 二.二项分布\n",
    "接下来我们对问题做一个的升级，假设我们一共抛了硬币$N$次，其中关于正面$x=1$出现次数$m$的概率分布称为二项分布，根据前面的推导，我们可以很容易写出它的分布：    \n",
    "\n",
    "$$\n",
    "Bin(m\\mid N,\\mu)=\\binom{N}{m}\\mu^m(1-\\mu)^{N-m}\n",
    "$$  \n",
    "\n",
    "其中，$\\binom{N}{m}=\\frac{N!}{(N-m)!m!}$，即$N$种组合中出现$m$次正面朝上的概率   \n",
    "\n",
    "#### 均值、方差\n",
    "如何按照定义求解会稍稍有些麻烦，由于二项分布可以看做$N$次独立的伯努利事件，那么**加和的均值等于均值的加和，加和的方差等于方差的加和**，我们知道$m=x_1+x_2+\\cdots+x_N$，所以：   \n",
    "\n",
    "$$\n",
    "E[m]=E[x_1]+E[x_2]+\\cdots+E[x_N]=N\\mu\\\\\n",
    "var[m]=var[x_1]+var[x_2]+\\cdots+var[x_N]=N\\mu(1-\\mu)\n",
    "$$\n",
    "\n",
    "#### 极大似然估计\n",
    "二项分布的极大似然估计求解其实和伯努利分布的求解一样，即在$N$次独立的伯努利实验中正面朝上出现的次数$m$所占的比例：    \n",
    "\n",
    "$$\n",
    "\\mu_{ML}=\\frac{m}{N}\n",
    "$$   \n",
    "\n",
    "那么，接下来让我们解决一下上面提到的过拟合问题，即为$\\mu$引入一个先验分布$p(\\mu)$，而对于二项分布比较有用的一个先验分布便是Beta分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三.Beta分布\n",
    "\n",
    "在构建先验分布时，我们往往想偷懒，因为如果先验分布比较复杂，再让它与似然函数相乘，那么后验的求解会更加困难，所以我们希望找到一种先验分布，它与似然函数相乘后也是易于分析的，最直接的一种方式就是构造一种与似然函数结构类似的先验分布，比如二项分布为某个因子与$\\mu^m(1-\\mu)^{N-m}$相乘的结果，那么我们同样构造一个类似结构先验分布，那么它们的乘积也会具有相同的结构，**即找到一个先验分布，让它与似然函数相乘后具有和先验分布相同的函数形式，这样的性质被称为共轭性**，那么我们可以假设我们的先验项中包含有$\\mu^a(1-\\mu)^b$这么一项，直观理解就是，我们先验分布中抛了$a+b$次硬币，其中$a$次正面朝上，$b$次反面朝上，由于需要做成一个分布，所以我们还需要找到一个归一化系数$\\beta$，使得：   \n",
    "\n",
    "$$\n",
    "\\int_0^1\\beta\\mu^a(1-\\mu)^bd\\mu=1\n",
    "$$   \n",
    "\n",
    "这里直接写出这个系数的表达式：   \n",
    "\n",
    "$$\n",
    "\\beta=\\frac{\\Gamma(a+b+2)}{\\Gamma(a+1)\\Gamma(b+1)}\n",
    "$$\n",
    "\n",
    "其中，$\\Gamma(\\cdot)$为Gamma函数，它的定义如下：   \n",
    "\n",
    "$$\n",
    "\\Gamma(x)=\\int_0^{\\infty}t^{x-1}e^{-t}dt\n",
    "$$  \n",
    "\n",
    "它看起来很奇怪，但它有一个很有用的性质：   \n",
    "\n",
    "$$\n",
    "\\Gamma(x)=(x-1)\\Gamma(x-1)\n",
    "$$  \n",
    "\n",
    "\n",
    "更多关于Gamma函数的由来，推荐博客：[《神奇的Gamma函数》](http://www.52nlp.cn/lda-math-%E7%A5%9E%E5%A5%87%E7%9A%84gamma%E5%87%BD%E6%95%B01)，接下来，我们就可以得到这个分布的形式了：   \n",
    "\n",
    "$$\n",
    "p(\\mu\\mid a,b)=\\frac{\\Gamma(a+b+2)}{\\Gamma(a+1)\\Gamma(b+1)}\\mu^a(1-\\mu)^b\n",
    "$$  \n",
    "\n",
    "但是呢...这个分布不是很美观，对超参做一个简单调整，令$a=a-1,b=b-1$，就得到了我们最终Beta分布的形式了：       \n",
    "\n",
    "$$\n",
    "Beta(\\mu\\mid a,b)=\\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)}\\mu^{a-1}(1-\\mu)^{b-1}\n",
    "$$  \n",
    "\n",
    "这里，$a,b$是先验分布中的超参数，关于不同$a,b$取值下的$\\mu$的分布如下图\n",
    "![avatar](./source/12_beta分布的超参数.png)\n",
    "可以发现一些简单的规律，比如若$a=1,b=1$，Beta分布则为$(0,1)$上的均匀分布，再比如随着$a,b$取值的增加，分布的“峰”会更尖，说明方差会越小，分布会越集中\n",
    "#### 均值、方差\n",
    "对于期望的求解，我们需要用到分布中自带的一个有用的性质：    \n",
    "\n",
    "$$\n",
    "\\int_0^1\\mu^{a-1}(1-\\mu)^{b-1}d\\mu=\\frac{\\Gamma(a)\\Gamma(b)}{\\Gamma(a+b)}\n",
    "$$   \n",
    "\n",
    "所以，均值：   \n",
    "\n",
    "$$\n",
    "E[\\mu]=\\int_0^1\\mu\\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)}\\mu^{a-1}(1-\\mu)^{b-1}d\\mu\\\\\n",
    "=\\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)}\\int_0^1\\mu^{(a+1)-1}(1-\\mu)^{b-1}d\\mu\\\\\n",
    "=\\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)}\\frac{\\Gamma(a+1)\\Gamma(b)}{\\Gamma(a+b+1)}\\\\\n",
    "=\\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)}\\frac{a\\Gamma(a)\\Gamma(b)}{(a+b)\\Gamma(a+b)}\\\\\n",
    "=\\frac{a}{a+b}\n",
    "$$   \n",
    "\n",
    "$E[\\mu^2]$的推导与$E[\\mu]$类似，只是上面积分项中$\\mu^{(a+1)-1}$变化为$\\mu^{(a+2)-1}$，所以： \n",
    "\n",
    "$$\n",
    "E[\\mu^2]=\\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)}\\frac{\\Gamma(a+2)\\Gamma(b)}{\\Gamma(a+b+2)}\\\\\n",
    "=\\frac{a(a+1)}{(a+b)(a+b+1)}\n",
    "$$   \n",
    "\n",
    "所以，方差：   \n",
    "\n",
    "$$\n",
    "var[\\mu]=E[\\mu^2]-(E[\\mu])^2=\\frac{ab}{(a+b)^2(a+b+1)}\n",
    "$$   \n",
    "\n",
    "#### 后验分布\n",
    "接下来看看后验概率分布长什么样子，我们将Beta分布和上面的二项分布相乘后，可以整理得到：    \n",
    "\n",
    "$$\n",
    "p(\\mu\\mid m,l,a,b)=\\frac{\\Gamma(m+a+l+b)}{\\Gamma(m+a)\\Gamma(l+b)}\\mu^{m+a-1}(1-\\mu)^{l+b-1}\n",
    "$$   \n",
    "\n",
    "这里，$N=m+l$，可以发现这也是一个Beta分布，即$Beta(\\mu\\mid a+m,b+l)$，显然这样的结果看起来也很make sense，我们可以通过它的均值来直观感受一下：   \n",
    "\n",
    "$$\n",
    "E[u]=\\frac{m+a}{m+a+l+b}\n",
    "$$\n",
    "\n",
    "#### 贝叶斯推断\n",
    "\n",
    "既然写出了后验分布的形式，那接下来就可以做预测了，比如下一次抛出正面的概率为：      \n",
    "\n",
    "$$\n",
    "\\int_0^1Bern(x=1\\mid\\mu)p(\\mu\\mid m,l,a,b)d\\mu=\\int_0^1\\mu p(\\mu\\mid m,l,a,b)d\\mu=E_{p(\\mu\\mid m,l,a,b)}[u]=\\frac{m+a}{m+a+l+b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 四.证明补充\n",
    "\n",
    "接下来对上面的内容做一些证明补充，第一个问题是若随机变量独立，那么它们的和的均值等于它们均值的和，以及它们的和的方差等于它们方差的和，第二个问题是Beta分布中归一化系数的推导，下面先看看第一个问题     \n",
    "\n",
    "（1）问题一，为了方便，这就可以推导两个独立的随机变量的情况，多个独立随机变量的推导类似，我们假设随机变量$x,y$独立，那么需要证明：    \n",
    "\n",
    "$$\n",
    "E[x+y]=E[x]+E[y]\\\\\n",
    "var[x+y]=var[x]+var[y]\n",
    "$$  \n",
    "\n",
    "下面证明一下，由于$x,y$独立，所以有$p(x,y)=p(x)p(y)$    \n",
    "\n",
    "所以，均值：    \n",
    "\n",
    "$$\n",
    "E[x+y]=\\int p(x,y)(x+y)dxdy\\\\\n",
    "=\\int p(x)p(y)(x+y)dxdy\\\\\n",
    "=\\int p(x)p(y)xdxdy+\\int p(x)p(y)ydxdy\\\\\n",
    "=[\\int p(x)xdx][\\int p(y)dy]+[\\int p(y)ydy][\\int p(x)dx]\\\\\n",
    "=\\int p(x)xdx+\\int p(y)ydy\\\\\n",
    "=E[x]+E[y]\n",
    "$$   \n",
    "\n",
    "方差：   \n",
    "\n",
    "$$\n",
    "var[x+y]=\\int\\int((x+y)-E[x+y])^2p(x,y)dxdy\\\\\n",
    "=\\int\\int((x-E[x])+(y-E[y]))^2p(x)p(y)dxdy（根据上面的推导，可以把E[x+y]拆开为E[x]+E[y]）\\\\\n",
    "=\\int\\int((x-E[x])^2+(y-E[y])^2+2(x-E[x])(y-E[y]))p(x)p(y)dxdy\\\\\n",
    "=[\\int(x-E[x])^2p(x)dx][\\int p(y)dy]+[\\int(y-E[y])^2p(y)dy][\\int p(x)dx]+2[\\int(x-E[x])p(x)dx][\\int(y-E[y])p(y)dy]\\\\\n",
    "=\\int(x-E[x])^2p(x)dx+\\int(y-E[y])^2p(y)dy+2(E[x]-E[x])(E[y]-E[y])\\\\\n",
    "=var[x]+var[y]\n",
    "$$   \n",
    "\n",
    "（2）问题二，Beta分布中归一化系数的推导其实就是验证下面的等式成立：    \n",
    "\n",
    "$$\n",
    "\\int_0^1\\mu^{a-1}(1-\\mu)^{b-1}d\\mu=\\frac{\\Gamma(a)\\Gamma(b)}{\\Gamma(a+b)}\n",
    "$$   \n",
    "\n",
    "以下推导[参考自>>>](https://blog.csdn.net/lanchunhui/article/details/75647076)，我们假设$t=x+y$，那么：   \n",
    "\n",
    "$$\n",
    "\\Gamma(a)\\Gamma(b)=\\int_0^\\infty e^{-x}x^{a-1}dx\\int_0^\\infty e^{-y}y^{b-1}dy\\\\\n",
    "=\\int_0^\\infty e^{-x}x^{a-1}[\\int_x^\\infty e^{x-t}(t-x)^{b-1}dt]dx（令y=t-x）\\\\\n",
    "=\\int_0^\\infty x^{a-1}[\\int_x^\\infty e^{-t}(t-x)^{b-1}dt]dx\\\\\n",
    "=\\int_0^\\infty e^{-t}[\\int_0^t x^{a-1}(t-x)^{b-1}dx]dt（交换积分顺序）\\\\\n",
    "=\\int_0^\\infty e^{-t}[\\int_0^1 (t\\mu)^{a-1}(t-t\\mu)^{b-1}td\\mu]dt（令x=t\\mu）\\\\\n",
    "=\\int_0^\\infty e^{-t}t^{a+b-1}dt\\int_0^1 \\mu^{a-1}(1-\\mu)^{b-1}d\\mu\\\\\n",
    "=\\Gamma(a+b)\\int_0^1 \\mu^{a-1}(1-\\mu)^{b-1}d\\mu\n",
    "$$  \n",
    "\n",
    "所以：   \n",
    "\n",
    "$$\n",
    "\\int_0^1\\mu^{a-1}(1-\\mu)^{b-1}d\\mu=\\frac{\\Gamma(a)\\Gamma(b)}{\\Gamma(a+b)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
